/* julia set fractals */

#include <stdio.h>
#include <stdlib.h>
#include "print_png.h"
#include <string.h>
#include <complex.h>
#include <omp.h>

// Maximum number of iterations
#define MAX_ITER 255
#define THRESHOLD 10
#define NUMTHREADS 32

#define POLY 1
#define SINE 2
#define COSINE 3

extern double wtime();

int main (int argc, char* argv[]){
  
  // Will be used to create the image
  char *values; 

  // Number of pixels in the generated square
  int size;
 
  // default values
  double xMin = -1.5;
  double xMax = 0.5;
  double yMin = -1.0;
  double yMax = 1.0;

  double t1, t2, t3, diff, print_time;

  char fileNameBase[50];
  char fileName[50];
  char line[20];

  double coef[4], temp; // corresponds to x^2 + x
  int i;

  // input file
  FILE * infile;

  // complex
  double real;// = -0.123;
  double imag;// = 0.745;
  complex double complexNum;

  if((argc!=4)){
    printf("Usage: julia.x $SIZE $OUTFILENAME $INFILE\n");
    exit(1);
  }

  infile = fopen(argv[3], "r");
  if(infile==NULL){
    printf("Error opening input file\n");
    exit(1);
  }
  else{
    printf("Reading from %s\n",argv[3]);
  }

  // complex point to calculate set for
  // fuck knows why i can't get fscanf to work right
  //fscanf(infile, "%f", &real);
  fgets(line, 20, infile);
  real = atof(line);
  fgets(line, 20, infile);
  //  fscanf(infile, "%f", &imag);
  imag = atof(line);
  complexNum = real+imag*_Complex_I;

  // box
  fgets(line, 20, infile);
  xMin = atof(line);
  fgets(line, 20, infile);
  xMax = atof(line);
  fgets(line, 20, infile);
  yMin = atof(line);
  fgets(line, 20, infile);
  yMax = atof(line);

  // polynomial
  for(i=3;i>=0;i--){
    fgets(line, 20, infile);
    //    fscanf(infile, "%f", &temp);
    coef[i]=atof(line);
  }

  fclose(infile);

  omp_set_num_threads(NUMTHREADS);
  
  if(1){
    printf("Read in: \n");
    printf("Point = %f + %fi\n", real,imag);
    printf("Box: xvals %f-%f, yvals %f-%f\n",xMin,xMax,yMin,yMax);
    printf("Poly: %fx^3 + %fx^2 + %fx + %f\n",coef[3],coef[2],coef[1],coef[0]);
  }

  size=atoi(argv[1]);
  
  strcpy(fileNameBase,argv[2]);

  values = (char*) malloc (sizeof(char)*size*size);
  if(values==NULL){
    fprintf(stderr,"Malloc failed\n");
    exit(1);
  }

  t1=wtime();
  
  calc_julia(xMin, xMax, yMin, yMax, values, size, complexNum, POLY, coef);
  strcpy(fileName, "POLY");
  strcat(fileName, fileNameBase);
  
  t2=wtime();
  
  print_png(values, size, size, fileName);
  /*calc_julia(xMin, xMax, yMin, yMax, values, size, complexNum, SINE, coef);
  strcpy(fileName, "SINE");
  strcat(fileName, fileNameBase);
  print_png(values, size, size, fileName);
  calc_julia(xMin, xMax, yMin, yMax, values, size, complexNum, COSINE, coef);
  strcpy(fileName, "COSINE");
  strcat(fileName, fileNameBase);
  print_png(values, size, size, fileName);*/
  
  //  t2=wtime();



  diff=t2-t1;

  printf("Time to run program = %f\n",diff);
  printf("%dx%d matrix\n",size,size);
  printf("xmin = %f, xmax = %f\n",xMin,xMax);
  printf("ymin = %f, ymax = %f\n",yMin,yMax);

  return 0;
}

/* Function calc_mb
   Calculates the number of iterations to determine whether or not a point
   is a part of the mandlebrot set
   Args xMin, xMax, yMin, yMax - map out the area for which the mandlebrot
     set is being calculated.  Not changed.
   Arg values - array in which the calculated values will be stored.
     Values are overwritten on output.
   Arg size - size of the nxn values array.  Not changed.
 */
int calc_julia(double xMin, double xMax, double yMin,
	       double yMax, char *values, int size,
	       complex double complexIn, int type, double* coef){

  int i, j, k;
  double x, y;
  //  double rn0, rn1, sn0, sn1;
  complex double point;
  complex double complexNum;
  
#pragma omp parallel for schedule(dynamic) private(j,k,x,y,complexNum,point)
  for(i=0; i<size; i++){
    for(j=0; j<size; j++){
      x = ((double)i/size)*(xMax-xMin) + xMin;
      y = ((double)j/size)*(yMax-yMin) + yMin;
      point = x + y * _Complex_I;
      complexNum = complexIn;
      for(k=0; k<MAX_ITER; k++){
     
	if(type==POLY){
	  point = coef[3] * cpow(point, 3.0) + coef[2] * cpow(point, 2.0) +
	    coef[1] * cpow(point, 1.0) + coef[0] + complexNum;
	  //	  complexNum = coef[3] * cpow(complexNum, 3.0) +
	  //  coef[2] * cpow(complexNum, 2.0) + coef[1] * complexNum + 
	  //  coef[0] + point;
	}
	if(type==SINE){
	  point=complexNum*csin(point);
	  // complexNum = point*csin(complexNum);
	}
	if(type==COSINE){
	  point = complexNum*ccos(point);
	  //complexNum = point*ccos(complexNum);
	}
	// rule of thumb cutoff that shows a point will not be in the set
	if(cabs(point)>THRESHOLD){
	  //	  printf("cabs = %f\n", cabs(complexNm));
	  break;
	}
	
      }
      values[j*size+i]=k;
    }
  }
}
